Plot standard errors (SE) from Local h2 estimation of Cross-Tissue and Tissue-Specific expression in GTex
my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
tislist <- scan(my.dir %&% 'GTEx_nine_tissues_spaces','c',sep="\n")
ct <- read.table(my.dir %&% 'GTEx.ranef.cross-tissue.h2_marginal.local_2015-03-23.txt',header=T)
nlist<-ct[1,2]
ts_allse <- ct %>% select(ensid,se)
colnames(ts_allse) <- c('ensid','CrossTissue')
for(i in 1:length(tislist)){
tis <- tislist[i]
tisdata <- read.table(my.dir %&% 'GTEx.resid.tissue-specific.h2_' %&% tis %&% '_marginal.local_2015-03-23.txt',header=T,sep="\t")
nlist <- c(nlist,tisdata[1,2])
tisse <- tisdata %>% select(ensid,se)
ts_allse <- inner_join(ts_allse,tisse,by='ensid')
tis <- gsub(' ','',tis) #replace spaces in tis string
tis <- gsub('(','',tis,fixed=T) #replace any parentheses in tis string
tis <- gsub(')','',tis,fixed=T)
tis <- gsub('-','',tis)
colnames(ts_allse)[2+i] <- tis
}
alllist<-c("Cross-Tissue",tislist)
samplesize<-data.frame(cbind(alllist,nlist))
colnames(samplesize)<-c("Tissue","n")
arrange(samplesize,desc(n))
## Tissue n
## 1 Cross-Tissue 450
## 2 Muscle - Skeletal 361
## 3 Whole Blood 339
## 4 Skin - Sun Exposed (Lower leg) 303
## 5 Adipose - Subcutaneous 298
## 6 Artery - Tibial 285
## 7 Lung 279
## 8 Thyroid 279
## 9 Nerve - Tibial 256
## 10 Heart - Left Ventricle 190
dim(ts_allse)
## [1] 17022 11
ts_a<-gather(ts_allse,"CrossTissue","Tissue",3:11)
colnames(ts_a) <- c('ensid','CrossTissue','TissueName','Tissue')
ggplot(ts_a,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Cross-Tissue SE') + xlab('Tissue-Specific SE') + ggtitle("GTEx Compare SE from Local h2 estimation")

Plot histograms of h2 and se Cross-Tissue and Tissue-Specific
histh2<-gather(ts_allh2,"ensid","Tissue",2:11)
colnames(histh2) <- c('ensid','Tissue','h2')
ggplot(histh2, aes(x = h2, fill = Tissue)) + geom_histogram() + facet_wrap(~Tissue,ncol=2) + ggtitle("GTEx Local h2") + guides(fill=FALSE)

ggplot(histh2, aes(x = h2, fill = Tissue)) + geom_histogram() + facet_wrap(~Tissue,ncol=2) + coord_cartesian(xlim=c(0.2,1),ylim=c(0,300)) + ggtitle("GTEx Local h2 > 0.20") + guides(fill=FALSE)

ggplot(histh2, aes(x = h2, fill = Tissue, color = Tissue)) + geom_density(alpha=0.3) + ggtitle("GTEx Local h2")

ggplot(histh2, aes(x = h2, fill = Tissue, color = Tissue)) + geom_density(alpha=0.3) + coord_cartesian(ylim=c(0,15)) + ggtitle("GTEx Local h2 ZOOM")

histse<-gather(ts_allse,"ensid","Tissue",2:11)
colnames(histse) <- c('ensid','Tissue','se')
ggplot(histse, aes(x = se, fill = Tissue, color = Tissue)) + geom_density(alpha=0.3)

se3<-inner_join(allse,ts_allse,by="ensid")
h23<-inner_join(allh2,ts_allh2,by="ensid")
seTW <- select(se3,1:11)
colnames(seTW)<-gsub('\\.x','',colnames(seTW))
gseTW <- gather(seTW,"Tissue","TissueWide",3:11)
dim(gseTW)
## [1] 153063 4
seTS <- select(se3,ensid,12:21)
colnames(seTS)<-gsub('\\.y','',colnames(seTS))
gseTS <- gather(seTS,"TissueName","TissueSpecific",3:11)
dim(gseTS)
## [1] 153063 4
gseALL <- cbind(gseTW,gseTS)
h2TW <- select(h23,1:11)
colnames(h2TW)<-gsub('\\.x','',colnames(h2TW))
gh2TW <- gather(h2TW,"Tissue","TissueWide",3:11)
dim(gh2TW)
## [1] 153063 4
h2TS <- select(h23,ensid,12:21)
colnames(h2TS)<-gsub('\\.y','',colnames(h2TS))
gh2TS <- gather(h2TS,"TissueName","TissueSpecific",3:11)
dim(gh2TS)
## [1] 153063 4
gh2ALL <- cbind(gh2TW,gh2TS)
ggplot(gh2ALL,aes(x=TissueSpecific,y=TissueWide)) +facet_wrap(~Tissue,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Tissue-Wide h2') + xlab('Tissue-Specific h2') + ggtitle("GTEx Local h2 estimates")

ggplot(gseALL,aes(x=TissueSpecific,y=TissueWide)) +facet_wrap(~Tissue,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Tissue-Wide SE') + xlab('Tissue-Specific SE') + ggtitle("GTEx Compare SE from Local h2 estimation")

ggplot(gh2ALL,aes(x=TissueSpecific+CrossTissue,y=TissueWide)) +facet_wrap(~Tissue,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Tissue-Wide h2') + xlab('Tissue-Specific + Cross-Tissue h2') + ggtitle("GTEx Local h2 estimates")

###make ordered point+CI h2 plots
gTW_h2 <- gather(h2TW,"Tissue.h2","TissueWide.h2",2:11)
gTW_se <- gather(seTW,"Tissue.se","TissueWide.se",2:11)
gTW_h2_se <- cbind(gTW_h2,gTW_se[,3])
colnames(gTW_h2_se) <- c("ensid","Tissue","h2","se")
ngenes<-nrow(gTW_h2_se)/length(unique(gTW_h2_se$Tissue))
gTW_h2_se <- gTW_h2_se %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) ) %>% arrange(Tissue,h2) %>% mutate(place=rep(1:ngenes,10),nonzeroCI=ymin>0)
ggplot(gTW_h2_se,aes(x=place,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + facet_wrap(~Tissue,ncol=2) + geom_pointrange(col='gray')+geom_point()+ggtitle('Tissue-Wide h2')+ theme(legend.justification=c(0,1), legend.position=c(0,1))+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))

gTS_h2 <- gather(h2TS,"Tissue.h2","TissueSpecific.h2",2:11)
gTS_se <- gather(seTS,"Tissue.se","TissueSpecific.se",2:11)
gTS_h2_se <- cbind(gTS_h2,gTS_se[,3])
colnames(gTS_h2_se) <- c("ensid","Tissue","h2","se")
ngenes<-nrow(gTS_h2_se)/length(unique(gTS_h2_se$Tissue))
gTS_h2_se <- gTS_h2_se %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) ) %>% arrange(Tissue,h2) %>% mutate(place=rep(1:ngenes,10),nonzeroCI=ymin>0)
ggplot(gTS_h2_se,aes(x=place,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + facet_wrap(~Tissue,ncol=2) + geom_pointrange(col='gray')+geom_point()+ggtitle('Tissue-Specific h2')+ theme(legend.justification=c(0,1), legend.position=c(0,1))+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))

##save data for later
write.table(h2TW,file=my.dir %&% "GTEx_Tissue-Wide_local_h2.txt",quote=F,row.names=F)
write.table(seTW,file= my.dir %&% "GTEx_Tissue-Wide_local_se.txt",quote=F,row.names=F)
write.table(h2TS,file=my.dir %&% "GTEx_Tissue-Specific_local_h2.txt",quote=F,row.names=F)
write.table(seTS,file=my.dir %&% "GTEx_Tissue-Specific_local_se.txt",quote=F,row.names=F)
TW<-inner_join(h2TW,seTW,by="ensid")
comboTW <- transmute(TW,EnsemblGeneID=ensid,h2.CrossTissue=CrossTissue.x,se.CrossTissue=CrossTissue.y,h2.AdiposeSubcutaneous=AdiposeSubcutaneous.x,se.AdiposeSubcutaneous=AdiposeSubcutaneous.y,h2.ArteryTibial=ArteryTibial.x,se.ArteryTibial=ArteryTibial.y,h2.HeartLeftVentricle=HeartLeftVentricle.x,se.HeartLeftVentricle=HeartLeftVentricle.y,h2.Lung=Lung.x,se.Lung=Lung.y,h2.MuscleSkeletal=MuscleSkeletal.x,se.MuscleSkeletal=MuscleSkeletal.y,h2.NerveTibial=NerveTibial.x,se.NerveTibial=NerveTibial.y,h2.SkinSunExposedLowerleg=SkinSunExposedLowerleg.x,se.SkinSunExposedLowerleg=SkinSunExposedLowerleg.y,h2.Thyroid=Thyroid.x,se.Thyroid=Thyroid.y,h2.WholeBlood=WholeBlood.x,se.WholeBlood=WholeBlood.y)
TS<-inner_join(h2TS,seTS,by="ensid")
comboTS <- transmute(TS,EnsemblGeneID=ensid,h2.CrossTissue=CrossTissue.x,se.CrossTissue=CrossTissue.y,h2.AdiposeSubcutaneous=AdiposeSubcutaneous.x,se.AdiposeSubcutaneous=AdiposeSubcutaneous.y,h2.ArteryTibial=ArteryTibial.x,se.ArteryTibial=ArteryTibial.y,h2.HeartLeftVentricle=HeartLeftVentricle.x,se.HeartLeftVentricle=HeartLeftVentricle.y,h2.Lung=Lung.x,se.Lung=Lung.y,h2.MuscleSkeletal=MuscleSkeletal.x,se.MuscleSkeletal=MuscleSkeletal.y,h2.NerveTibial=NerveTibial.x,se.NerveTibial=NerveTibial.y,h2.SkinSunExposedLowerleg=SkinSunExposedLowerleg.x,se.SkinSunExposedLowerleg=SkinSunExposedLowerleg.y,h2.Thyroid=Thyroid.x,se.Thyroid=Thyroid.y,h2.WholeBlood=WholeBlood.x,se.WholeBlood=WholeBlood.y)
write.table(comboTW,file= my.dir %&% "GTEx_Tissue-Wide_local_h2_se.txt",quote=F,row.names=F,sep="\t")
write.table(comboTS,file=my.dir %&% "GTEx_Tissue-Specific_local_h2_se.txt",quote=F,row.names=F,sep="\t")